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A fluid droplet may exhibit self-propelled motion by modifying the wetting properties of the 
substrate. We propose a novel model for droplet propagation upon a terraced landscape of ordered 
layers formed as a result of surface freezing driven by the contact angle dependence on the terrace 
thickness. Simultaneous melting or freezing of the terrace edge results in a joint droplet-terrace 
motion. The model is tested numerically and compared to experimental observations on long-chain 
alkane system in the vicinity of the surface melting point. 
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Motion of mesoscopic liquid droplets is a challenging 
problem both in view of numerous technological applica- 
tion in surface treatment, microfluidics, etc., and fun- 
damental questions arising on the borderline between 
macroscopic and molecular physics. Different scenarios 
of droplet motion are determined by liquid-substrate in- 
teractions, and may incorporate surface phase transitions 
and chemical reactions, as well as more subtle modifica- 
tion of physical properties in interfacial regions. One 
can distinguish between three classes of behavior: pas- 
sive, interacting, and active. A passive droplet gains 
mobility due to imposed forces, e.g. temperature gra- 
dients Q or substrate heterogeneity 0]. Motion of in- 
teracting droplets is mediated by fluxes through a thin 
precursor layer 0, Q ■ Finally, active droplets may pro- 
pel themselves by modifying the substrate either through 
surfactant deposition at the three-phase contact line |5j 
or through chemical reaction proceeding directly on the 
substrate at the foot of the droplet 0, Q • 

A new type of self-propelled motion discovered recently 
in experiments with long-chain alkanes (C n H2 n +2) is 
associated with surface phase transitions creating a ter- 
raced immobilized layer between the fluid and substrate, 
as shown schematically in Fig. ^ The system includes 
(a) a disordered bulk liquid alkane droplet; (b) one or 
more ordered (smectic A) alkane layers formed as a re- 
sult of surface freezing; (c) a molecularly thin disordered 
precursor layer. The thickness ratio l/d^> 1 of the smec- 
tic and precursor layers is determined by the aspect ratio 
of the alkane molecule. The plateau height is H = Nl, 
where ./V > 1 is an integer. A similar situation may arise 
in layered adsorption, leading to the formation of ordered 
immobilized molecular layers also in the case of the as- 
pect ratio l/d ~ 0(1). 

Due to a difference in molecular interaction strengths 
between the bulk fluid and the smectic and the sub- 




FIG. 1: A scheme of a droplet on a terraced smectic layer. 
I and d denote, respectively, the molecular dimensions along 
and across the long molecular axis; H — Nl is the terrace 
height, N is an integer that denotes number of layers, h is 



the bulk droplet height, and 6 
upper and lower plateaus. 



are the contact angles on the 
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strate, the contact angle of the bulk droplet depends on 
the number of smectic layers, and therefore the droplet 
is expected to move when placed on a terrace as in 
Fig. ^ Moreover, as temperature is varied, the terrace 
surface freezing process may proceed in two ways, de- 
pending on whether it is limited by material supply or 
removal of latent heat. The first mechanism involves slow 
spreading, with the smectic layer growing sidewise, be- 
ing augmented by fluid molecules migrating from the bulk 
droplet through the precursor [t| . The second mechanism 
is fast, and involves terrace growth synchronous with the 
droplet motion Melting, being unconstrained by ma- 
terial supply, proceeds by the second mechanism in re- 
verse. In this Communication, we suggest a model of self 
propelled droplet motion accompanied by surface freez- 
ing or melting on terraced landscape. 

We adopt the lubrication approximation, which ac- 
counts for different scalings in the vertical and the hor- 
izontal directions |10|. The approximation is applicable 
in a liquid film with a large aspect ratio, when the in- 
terface is weakly inclined and curved. The scaling is 
consistent if one assumes d z ~ 0(1), V ~ O(e) <C 1, 
dt ~ C(e 2 ), where V is the two-dimensional gradient in 
the horizontal plane. This scaling also implies a small 
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contact angle, 8 ~ 0(1) and results in a different or- 
der of magnitude of the vertical and horizontal veloci- 



0(e). As a consequence, 



ties, v z ~ 0(e ), ~ v, 
the pressure or, more generally, a driving potential W is 
constant to 0(e 2 ) across the layer in z direction. The 
governing equation for the droplet height h following the 
mass conservation condition reads: 

d t h = -W-j, j = rj- 1 k(h)VW, W = Il-ae 2 V 2 h, 

(1) 

where j is the mass flux, 77 is the dynamic viscosity, k(h) is 
the mobility coefficient, tr is the surface tension, and II is 
the disjoining potential due to interaction with the solid 
support (including both substrate and smectic layers). 

Computation of II is the key component of the model. 
We assume that all interactions are of the van der Waals 
type with the hard core potential V(r) = 00 at r < d, 
V(r) = —Ajr~ 6 at r > d, and differ by interaction con- 
stants Aj only. Since the motion is, on the one hand, 
caused by the difference in contact angles on the two 
sides and, on the other side, is driven by excess free en- 
ergy of either freezing or melting, the difference between 
the liquid-terrace (At) and liquid-substrate (A s ) interac- 
tion constants should change sign when the temperature 
passes the surface melting point. Equilibrium contact an- 
gles can be expressed through the interaction constants 
by integrating Eq. Q for an infinite bulk fluid in equi- 
librium as explained below. 

For a fluid on top of a flat homogeneous plateau, z > H 
(see Fig.^), the free energy per unit area can be written, 
in a local density functional approximation |l l| . as 

7 = / n(z) lf(n)- / Q(*-C)n(*)dC+ 
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Here n(z), n t , n s are the fluid, terrace, substrate par- 
ticle densities and f(n) is free energy per particle of a 
homogeneous fluid. The first term in the integrand is 
the free energy per particle in the homogeneous fluid; 
the second term compensates lost fluid-fluid interactions 
in the domain z < H which are included in f(n); the 
third term accounts for the inhomogeneous part of fluid- 
fluid interactions; the last two terms represent the fluid- 
terrace and fluid-substrate interactions. All interactions 
are described by the same hard core interaction poten- 
tial differing only by interaction strength, A (fluid- fluid) , 
A t = at A (fluid-terrace) and A s = a s A (fluid-substrate). 
The interaction kernel Q(C) lumping intermolecular in- 
teraction between the layers z — constant [ll| is ex- 
pressed then as Q(() — ^nA(~ 4 at C > d. 

Since the precursor layer is of molecular thickness, 
the chemical potential shift is computed differently in 
the bulk and precursor regions; this is unlike other self- 



propelled active drop models 7) where a macroscopic pre- 
cursor layer was presumed. In the bulk region z > H + d 
the chemical potential shift fi(h) — no from the equilib- 
rium value in the bulk fluid, fj, = /xo, depends on the 
fluid thickness h and coincides with the disjoining po- 
tential, H(h) = dhj [HI]- Neglecting the vapor density, 
as well as density variation in a molecularly-thin inter- 
facial layer, we can apply the sharp interface approxi- 
mation [13, assuming the fluid density to be constant, 
n = rio at H + d < z < h, where n$ is the equilibrium 
fluid particle density at /1 = /io, n = at z > h. Defining 
j(h) by Eq. J5| with the upper integration limit over z 
replaced by h and the homogeneous part excluded, we 
compute 
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where \ = a t n t /n - 1 and \a = (a s n s - a t n t )/n are 
dimensionless Hamaker constants for fluid-terrace and 
terrace-substrate interfaces. 

The precursor film is assumed to be of a constant 
molecular thickness d, but the liquid density is allowed 
to vary there, and is determined by minimizing the grand 
ensemble thermodynamic potential T — 7 — /ijndz. 
The disjoining potential is identified here with the shift of 
chemical potential per unit volume II (n) = n[fi(n) — /io] 
relative to the equilibrium value fio as a function of the 
local value of n (shifted from its bulk equilibrium value 
under the action of the terrace and substrate). It is de- 
termined by the Euler-Lagrange equation derived from 
the integrand of J5J for z = H + d: 
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The mobility coefRcient, k(h), is also computed sepa- 
rately in the bulk and precursor regions and matched at 
the precursor thickness. In the bulk region, Stokes flow 
with a kinematic slip condition [l4j | is assumed, while 
in the precursor domain the mass transport is presumed 
to be purely diffusional. This yields the mobility coeffi- 
cient 



k(h) = i x2 ( h -H) + ^[ h -( H + d )l at h > H + d; 
[ ' \ X 2 d at h < H + d , 

(5) 

where A = J Dn/n^kBT ~ O(d) is the effective slip 
length; D is surface diffusivity, fcg is Boltzmann con- 
stant, and T is temperature. 

The motion of a droplet placed on terraced landscape 
as in Fig. ^ can be attributed to a difference in equilib- 
rium contact angles at the upper (H + ) and lower (H~) 
terraces. The rescaled angles can be calculated for \ < 0, 
\x\ <§C 1 by integrating the static equation W = [TT| . 
which reduces to ae 2 h xx = II. In the limit h — > 00 we 
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obtain 
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(6) 

where ho « i? ± + <i. The direction of the droplet motion 
is determined solely by the effective terrace-substrate in- 
teraction, i.e. by the sign of Xa- the droplet either as- 
cends for Xa > or descends for Xa < until equi- 
librium is reached. The equilibrium condition 9 + = 6~ 
is satisfied either by H + = H~ or Xa = 0. The for- 
mal small parameter of the lubrication approximation 
can b e defined by se tting 9 = 1 for Xa = 0, which yields 
e ~ \/|x] An 2 /(ad 2 ). Since a ~ n 2 A/d 2 , a good estimate 

is e - VTxT- 

A decrement of contact angles should be preserved to 
maintain droplet propagation. This is possible when the 
terrace edge is also allowed to move. The terrace motion 
due to surface freezing or melting was observed in the 
experiment when ambient temperature T was varied 
in the vicinity of the surface freezing point T m . When 
the terrace is at the foot of a liquid droplet, the melting 
or freezing rate is limited by the heat flux q required to 
supply or remove the latent heat L, so that Lpv = g« 
K(T — T m )/(h — H + ), where /C is thermal conductivity 
and p is density (assumed to be equal for both liquid and 
the frozen terrace layer). The approximate expression for 
the heat flux (directed almost normally to the substrate 
or terrace) corresponds to the lubrication approximation. 
The form of this relation defining the edge position x is 

_ da; _ JC(T - T m ) 
V ~ dt ~ Lp (h - H+) ' 

To reproduce joint droplet-terrace dynamics observed 
in H, we have carried out dimensionless ID numerical 
computations of Eqs. 0J, JJJ). The new dimensionless 
variable forms are: h — h/d, £ = xe/d, t = te a/(drj) 
and II = nd/(e 2 cr). The particle densities are scaled by 
1/6, where b — 27rd 3 /3 is the excluded volume so that 
the respective dimensionless equations are: 
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FIG. 2: Numerical solution of the fluid and the terrace ac- 
cording to Eq. ©, showing the melting process (a-d) at re- 
spective time steps (from top to bottom: 50, 210, 400, and 
650). The horizontal range is £ = [0, 120] and the vertical 
range is h = [0, 20]. The dotted line marks the droplet height 
above the terrace edge, h c . The inset shows the dependence 
of h c — H + as a function of time and its relaxation to an equi- 
librium value h c — H + ~ 4. Parameters: x = — 0.3, /3 = 15, 
A = Vs, I — 2, H + — 21, H~ = l, Xa = -10 and A = -4. 



(3 = A/(k B Td 6 ) and jtt = n b/{k B T). The same no- 
tation is retained for the dimensionless variables, and 
Eq. 10 remains without change, except replacing d — > 1. 
The density in the precursor domain, n is transformed 
to effective height as n = rio(h — H). We adopted the 
explicit spectral method and by doubling the grid size 
impose reflecting boundary conditions. The initial state 
in each computation includes a droplet with its maxi- 
mum placed above the terrace edge and a precursor film 
of unit thickness (see Fig. QJ. The parameter A defining 
the ratio of characteristic velocities of the edge to droplet 
motion is of 0(1) when the temperature difference is in 
the range of O(10 _3 )[°K]- 

Synchronous droplet-terrace motion under melting 
conditions is shown in Fig. [3 This joint propagation 
can be explained in terms of terminal velocity of the ter- 
race. While the terrace is below the droplet, the droplet 
velocity is determined by the difference in contact angle 
values, according to ©■ On the other hand, the terrace 
velocity (for a fixed A) depends solely on the thickness 
of the liquid layer above the terrace edge. At the start, 
the droplet moves to the left, while the terrace remains 
almost stationary because of slow transport through a 
thick layer, as shown Fig. |2a-b). As the fluid height 
above the edge decreases, the terrace gains speed [see 
Fig. |2£b-c)] , until it reaches an "equilibrium" position, 
such that the point at the droplet interface just above the 
edge where the thickness is h c moves with the same speed 
v = A/(h c — H + ). The stable position should lie near 
the trailing edge; then, if the terrace moves faster than 
the droplet, the liquid layer thickness above it increases 
and the terrace decelerates. As a result, the value h = h c 
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FIG. 3: Numerical solution of Eq. JHj, showing the freezing 
process (a-e) at respective time steps (from top to bottom: 0, 
410, 581.85, 581.88 and 581.92). The horizontal range is £ = 
[0, 120] and the vertical range is h — [0, 12]. The dotted line 
in (c-e) marks the droplet position according to its maximal 
height. The inset shows the dependence of h c — H + on time 
in the vicinity of the critical droplet volume. Parameters: 
X = -0.3, (3 = 15, A = y/3, I = 2, H+ = 21, H~ =l, Xa = 1.5 
and A = 4. 

remains constant, as seen in Fig. Hfc-d)] and more pre- 
cisely in inset of Fig.[2d). This dynamic feedback allows 
the droplet and the terrace to synchronize their motion. 
As implied by © , we found that the synchronous propa- 
gation velocity depends on the layers thickness H + = Nl 



and will be discussed elsewhere |15J. 

In the freezing case, the droplet volume decreases, since 
the total mass of the system is conserved and the fluid is 
solidified. The droplet and the terrace may still jointly 
propagate as long as the droplet height is relatively large 
compared to the precursor thickness, as shown in Fig^a- 
c) . In a such motion the droplet and the terrace preserve 
the equilibrium height h = h c [see Fig. E^b-c)]. As the 
droplet volume decreases below the equilibrium height 
h c , the terrace propagates faster than the droplet and 
runs out to its leading edge [see Fig. EJd-e)]. Following 
this, the motion stops, since further terrace propagation 
is limited by slow material supply through the precursor, 
and the droplet is left in an equilibrium state on the top 
of a flat smectic layer 8] . This behavior is also presented 
in the inset of Fig. 13 As the terrace passes the maximum 
droplet 's height, the critical value h c — H + decreases to 
unity (i.e. the precursor thickness). The droplet velocity 
at the same time drops to zero, while the terrace velocity 
(dashed line) jumps abruptly. Similar behavior has been 
also observed experimentally [jj. 

We have proposed a model for self-propelled droplets 
on top of a terraced landscape driven by surface freezing 
or melting. The numerical estimates show the charac- 
teristic terrace edge velocity v ~ O(10 2 ) {[xm/sec] close 
to the experimental data |8( at temperature variations 
around the surface melting temperature \T — T m \ ~ 
O(10- 3 )[°K] and h c ~ 0(d) ~ 0.1 [ram]. 
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